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Macroscopic yielding in jammed solids is accompanied by a non-equilibrium 
first-order transition in particle trajectories 
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We use computer simulations to analyse the yielding transition during large-amplitude oscillatory 
shear of a simple model for soft jammed solids. Simultaneous analysis of global mechanical response 
and particle-scale motion demonstrates that macroscopic yielding, revealed by a smooth crossover in 
mechanical properties, is accompanied by a sudden change in the particle dynamics, which evolves 
from non-diffusive motion to irreversible diffusion as the amplitude of the shear is increased. We 
provide numerical evidence that this sharp change corresponds to a non-equilibrium first-order 
dynamic phase transition, thus establishing the existence of a well-defined microscopic dynamic 
signature of the yielding transition in amorphous materials in oscillatory shear. 
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I. INTRODUCTION 

A major effort in soft condensed matter physics con¬ 
cerns the design of materials with well-controlled me¬ 
chanical properties El- Rheology thus represents a cen¬ 
tral probe and oscillatory shear measurements at finite 
frequency u> are among the most commonly performed 
mechanical tests (2j. In this approach, a harmonic de¬ 
formation is applied and the stress response measured, 
or vice-versa. In the linear regime, the complex shear 
modulus G*(w) = G'(w) + i G"{u) provides information 
about the nature and strength of the material at a given 
frequency Q, while microscopic relaxation processes can 
be probed by varying the frequency. At larger amplitude, 
non-linear mechanical properties are accessed. 

This approach is well-suited for amorphous materials, 
which often display non-trivial response spectra in the 
linear regime where they behave as soft elastic solids, 
but flow as the amplitude of the forcing is increased be¬ 
yond a “yielding” limit Q-Q- Whereas the storage mod¬ 
ulus G'(w) dominates the elastic response at small defor¬ 
mation amplitude, irreversible plastic deformations oc¬ 
cur post yielding where the loss modulus G"(w) instead 
dominates. Oscillatory shear experiments have been per¬ 
formed in a wide range of soft condensed matter systems 
across yielding, such as granular particles EH0 , emul¬ 
sions fl3 - fl5j |. colloidal suspensions ITgI-ITbI and gels 0, 
as well as in computer simulations [20l424| . 

In experiments, the change from elastic to plastic re¬ 
sponse in macroscopic mechanical properties is often de¬ 
scribed as a “yielding transition”, even though yielding 
appears as a smooth crossover whose location cannot 
be unambiguously defined [25]. Interestingly, recent ex¬ 
periments have provided evidence that this macroscopic 
crossover corresponds to a qualitative change in particle 
trajectories @, [lj EM, EM EMUM- As expected physi¬ 
cally, particles are essentially arrested in the undeformed 
solid, but can diffuse due to irreversible plastic rearrange¬ 
ments occurring at larger amplitude. There is, how¬ 


ever, no consensus about the nature of this crossover, 
which has been described either as a smooth change 0, 
as a relatively sharp crossover [121. or as a continuous 
non-equilibrium phase transition [18(. The latter con¬ 
clusion builds a qualitative analogy with the continu¬ 
ous irreversibility transition observed in low-density sus¬ 
pensions 10 . which has been actively studied in 

computer simulations [23j, [0 , and attempts to borrow 
conce pts from the field of non-equilibrium phase transi¬ 
tions [30|j32|. In addition recent experiments argue that 
yielding cor resp onds to a change in the microstructure of 
the system [27J , by opposition to the dynamic properties 
discussed here. A clear connection between these micro¬ 
scopic changes and the macroscopic rheology is lacking. 

Here, we use a model of a jammed material composed 
of non-Brownian repulsive spheres to investigate the na¬ 
ture of the yielding transition at the particle-scale level, 
in the simple situation where thermal fluctuations and 
hydrodynamic forces play no role. We reproduce stan¬ 
dard mechanical signatures of macroscopic yielding un¬ 
der oscillatory shear and obtain two key results regarding 
particle trajectories. First, we show that the onset of par¬ 
ticle diffusion in steady state is extremely sharp and oc¬ 
curs at a well-defined shear amplitude, which unambigu¬ 
ously locates the yielding transition. Second, we find that 
particle diffusivity emerges discontinuously at yielding, 
thus demonstrating that yielding corresponds to a non¬ 
equilibrium first-order phase transition. These findings 
differ qualitatively from earlier suggestions of a continu¬ 
ous irreversibility transition EM EH , but seem to agree 
very well with recent experimental findings 00. We 
also show that this transition is dynamic in nature, and 
is not accompanied by discontinuous structural changes. 


II. MODEL AND NUMERICAL TECHNIQUES 

We consider soft repulsive non-Brownian particles in a 
simple shear flow geometry. We perform standard over- 
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damped Langevin dynamics simulations of a well-known 
model of harmonic particles in three dimensions [33i . f34j | , 
using an equimolar binary mixture of small and large par¬ 
ticles with diameter ratio 1.4. The equations of motion 
read 


-7 {t)Vi(t)e x 


£ 


dU(r t 
dr a 


= 0 , 


(1) 


where £ s is a friction coefficient, fij = (xij, yij , Zij) = 
(xj—Xi,yj — yi,Zj — Zi), e x = (1,0,0) and j(t) is the shear 
rate. For particles i and j having diameters and a,j , the 
pair potential reads [/(r^) = | (1 — r^/ay) 0 (aij — ry,-), 
where e is an energy scale, a,;, = (cq + a j)/ 2 , a denotes 
the diameter of small particles, and 0(x) is the Heaviside 
function. The unit length is a, the unit time tq = a 2 £ s /e, 
the unit energy e, and so the unit stress is e/a 3 . 

We apply a harmonic deformation using Lees-Edwards 
periodic boundary conditions (34j |, the strain evolving as 
7 (f) = 7 o[l — cos(ujt)], where 70 is the amplitude of the 
imposed shear strain and u> = t/T is the frequency of 
the oscillation. The period T is chosen very large to 
be close to a quasi-static protocol; we use T = 10 4 to- 
We have checked that our results do not qualitatively 
depend on this choice. We work at constant packing 
fraction <p = 0.80, much above the jamming density 
ipj ~ 0.647 & We checked that our results are repre¬ 
sentative of the entire jammed phase, p > <pj, but yield¬ 
ing could be more complicated in the limit <p —>• tpj where 
the system looses rigidity [l2j]. The different regime 
ip < ipj, where a yield stress does not exist, was anal¬ 
ysed befor e |29fl . We solve Eq. (JT]) with a modified Euler 
algorithm [.'1 11 using a discretization timestep At = 0.1 t 0 . 
Numerical stability and accuracy were carefully checked. 
To investigate finite-size effects, we perform simulations 
with four different sizes N = 300, 1000, 3000, and 10000, 
where N is the total number of particles. All simulations 
start from fully random configurations. We analyse both 
the transient regime after shear is started, and steady- 
state measurements. To improve the statistics, we per¬ 
form at least 4 independent runs for each pair ( 70 , N). 

At the macroscopic level, our main observable is the 
time-dependent response of the xy-component of the 
shear stress, defined by the usual Irving-Kirkwood for¬ 
mula [34jj: cr(t) = —y^2 x ijFij, where V is the volume 
and F-j represents the y-component of the force Fij . The 
kinetic part of the stress is fully negligible in the present 
situation of low-frequency oscillatory shear. To analyse 
the rheological response in steady state, we fit the time 
series of the shear stress to a sinusoidal form, 


a(t) = — (To cos (wt + 5), ( 2 ) 



FIG. 1. (a) Evolution of storage and loss moduli, Eq. 0, 
with strain amplitude for different system sizes. The crossing 
of G' and G" at (dashed lines) defines a characteristic 
strain scale, (b) Evolution of amplitude and phase of the 
stress, Eq. 0, with strain amplitude. The amplitude has a 
maximum at 7 P i (dashed lines), (c) System size dependence 
of all characteristic strain amplitudes. 


the more conventional quantities G'(u>) and G"{u) using 

G'{lo) + \G"{pj) = —e iS . (3) 

7o 


III. SMOOTH CROSSOVER IN MACROSCOPIC 
RHEOLOGY 


where <to is the amplitude of the first harmonics at fre¬ 
quency w, and S is the phase difference between strain 
and stress. In practice, a Q and 5 are obtained by fitting 
Eq. @ to steady state data lasting about lOOT. Alterna¬ 
tively, we can transform the two parameters (<toj< 5) into 


In Fig. [lja), we show the evolution with the strain 
amplitude 70 of the storage and loss moduli at fixed fre¬ 
quency uj measured in steady state. At very low 70 , G' 
dominates the response, G'/G" ~ 10, indicating that 
the system responds in the linear regime as a soft elastic 
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solid. As 70 is increased, the moduli first evolve slowly 
for 7 q < 0 . 1 , where little plastic rearrangements are pro¬ 
duced. As 70 increases further, we observe a crossing of 
G' and G" at 7 X « 0.15 (dashed lines), so that dissipa¬ 
tion dominates 70 > 7 x • These mechanical properties 
reproduce well-known behaviour 0,1 and validate our 
numerical approach. We notice further that they dis¬ 
play virtually no finite-size effects. In this representa¬ 
tion, 7 X ss 0.15 appears as the most relevant strain scale 
to characterize yielding, although a smooth crossover can 
be qualitatively detected near 70 « 0 . 1 , where the 70 - 
dependence of the moduli becomes somewhat steeper. 

In Fig. mb) we replot this rheological evolution in the 
alternative representation offered by (oo,^)- At small 
70 , the phase <5 is very small while the stress amplitude 
increases linearly, cto ~ G'y 0 , as expected for reversible 
elastic deformations in a solid. Near 70 w 0.1 two changes 
are observed. First, cr 0 ceases to be linear and displays 
an overshoot when 70 is increased, signaling that plastic 
events take place. We define the onset of plastic events, 
7 P i, as the location of the stress overshoot [dashed lines in 
Fig.[ljb)]. Our interpretation for 7 p i is reinforced by the 
evolution of the phase <5 in Fig.QJb), which grows steadily 
above 7 P i, indicating the onset of dissipation. Note that 
the crossing of G' and G" at y x has no obvious relevance 
in this representation where it simply corresponds, by 
definition, to the strain scale where <5 = 7 r/ 4 . In Fig.[ljc) 
we confirm that 7 P i and y x display virtually no system 
size dependence, but that they differ quantitatively. 

Whereas y x is frequently quoted as “the” yielding 
point in the literature, a stress overshoot also serves to 
identify yielding in shear-start experiments [25 . A stress 
overshoot is reported in some oscillatory shear experi¬ 
ments [li|, but is absent in others [TH, [l3j] • A possible 
explanation is that experiments are typically performed 
at somewhat larger frequencies, where additional contri¬ 
butions to the shear stress (lubrication forces, hydrody¬ 
namic effects) might hide this behaviour. In addition, we 
have confirmed that the overshoot disappears and is re¬ 
placed by a monotonic increase seen in experiments [l3| 
when we use a substantially larger frequency, typically 
uj > 10 -3 . We emphasize that the presence of the dy¬ 
namic transition discussed below is independent of the 
existence of the stress overshoot reported in Fig. [TJb). 


IV. SHARP TRANSITION IN MICROSCOPIC 
DYNAMICS 
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FIG. 2. (a) Transient evolution of particle displacement 

A r(t,T), Eq. (f4~fl . for various strain amplitudes for N = 3000. 
For 70 < 0.090, the displacements drop to zero after a 
timescale ra. (b) The time averaged displacement for one 
cycle (Ar(T)) for various strain amplitudes and system sizes. 
The stable and the metastable data are plotted with filled 
and open symbols respectively. The vertical dashed lines in¬ 
dicate 7 c for each system size, (c) The average ra diverges 
algebraically, ra ~ ( 7 a yn — 70 ) interpreted as the diverg¬ 
ing lifetime of a metastable state near a discontinuous phase 
transition. The vertical dashed lines indicate 7 a yn for each 
system size. 


We now turn to the evolution of single particle dynam¬ 
ics. A first natural dynamic observable is the averaged 
particle displacement after one deformation cycle |l2U24j . 

Ar(t,T) = ^^|r i (i + T)-r i (t)|, (4) 

j 

where t is the time since shear is applied. In Fig. [21(a), 
we show how Ar(f, T) evolves in the transient regime 


for various amplitudes of the applied deformation. For 
small 70 , particle displacements decay rapidly to zero. In 
the elastic solid at small amplitude, there are rare rear¬ 
rangements taking place before the system settles near a 
stable energy minimum where particles have nearly peri¬ 
odic motion (or quasi-periodic motion with a period that 
is a multiple of T, as reported before (H, l3fij|). As 70 
is increased, it takes more and more time for Ar(f, T ) 
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to eventually vanish. When 70 is larger than 7 ~ 0.095, 
the average particle displacement never vanishes in the 
explored time window, but instead fluctuates around a 
well-defined finite value, which increases with 70 . This 
regime corresponds to irreversible, non-periodic particle 
trajectories. In Fig. [2jb), we plot the time averaged dis¬ 
placement for one cycle (A r(T)) in steady states for vari¬ 
ous strain amplitudes and system sizes. From (Ar(T)), a 
very clear discontinuous jump is observed between the ir¬ 
reversible and reversible states near 7 C . Very close to the 
transition, the displacements exhibit fluctuations around 
a well-defined value both above and below 7 C . Whereas 
these fluctuations are infinitely long-lived above j c , they 
are only metastable below 7 C , before the system finds 
a reversible state where the displacements become very 
small. We report the value of (A r(T)) for N = 10000 
for these metastable states in Fig. [2f b). Overall, these 
fluctuations appear qualitatively distinct from the alge¬ 
braic decay observed close to continuous irreversibility 
transitions Q, and are much closer to the phenomenol¬ 
ogy observed near discontinuous, first-order phase transi¬ 
tions where metastable phases can be observed over long 
times. In particular, it appears impossible to describe 
the decrease of (A r(T)) with a continuous vanishing at 
the critical value of 7 d yn - 

Stronger evidence of such a phase transition is ob¬ 
tained from the evolution of the average lifetime of the 
metastable irreversible phase Td, as depicted in Fig. [31c) 
for various system sizes. These data confirm that Td in¬ 
creases rapidly close to 70 ~ 0.1. For larger 70 , trajec¬ 
tories remain irreversible. By contrast to the rheology, 
a clear system size dependence is observed, larger sys¬ 
tems take more time to settle in a global energy mini¬ 
mum. A diverging lifetime is typically observed close to 
non-equilibrium phase transitions (§1, [301, [3l|, and was re¬ 
ported before [H, [HJ . Such divergence is expected for 
both a continuous or discontinuous transition (see [37ll38j ] 
for recent examples). Solid lines in Fig. 0)c) represent an 
emp irical fit, Td ~ ( 7 dyn — 7 o) _a , with a ~ 2.1 — 3.0 
[39(, suggestive of a divergence of Td when approaching 
the dynamic transition at 7 d yn - We notice that finite size 
effects can be felt even when the system is not very close 
to the critical value, such that we cannot observe the ex¬ 
pected saturation of Td to a finite value as N —> 00 . The 
system size dependence of 7 d yn reported in Fig. Oja) is 
very modest and seems to extrapolate to a finite value, 
7 dyn ~ 0.095, as N —> 00 . 

We now characterize the steady-state irreversible dy¬ 
namics at large 70 using the mean-squared displacement, 

1 N 

i=i 

where the brackets indicate a time average. The results 
are displayed in Fig. [31(a) for N = 10 4 , for time delays 
commensurate with the period. The dynamics is diffusive 
at long times, (A r(t) 2 ) ~ 6 Dt 1 where D is the diffusion 
constant. We represent D (in units of a 2 /T) for various 



FIG. 3. (a) Mean-squared displacements, Eq. ([5]). for N = 10 4 
and various strain amplitude become diffusive at long time 
when 70 > 7 c ~ 0.085. (b) The diffusion constant D decreases 
modestly as 70 decreases towards a critical value -y c where 
it drops discontinuously to zero. The vertical dashed lines 
indicate 7 C for each system size and the black one represents 
the large-A extrapolation of y c « 0.0885 performed in the 
inset. 


70 and N in Fig.[3Ib). As expected, D = 0 below a criti¬ 
cal value 7 c, corresponding to the phase characterized by 
quasi-periodic particle trajectories, and it increases with 
70 above j c . Both D and A r(t,T) in Eq. (J4J could serve 
as order parameters for the transition. By measuring 
7 C for various system sizes, we observe a modest change 
with system size, see inset of Fig. [31b), suggestive of a 
finite limit y c ss 0.0885 for N —> 00 . The functional form 
of our extrapolation should be confirmed by additional 
larger scale simulations. 

A striking finding in Fig. [31b) is the finite amplitude of 
the diffusion constant at the transition. Near continuous 
irreversibility transitions, D decreases by several orders 
of magnitude and scales algebraically as j c is approached 
from above [3, SOj] . We observe instead a modest decrease 
of D, followed by a sudden jump to zero, which is robust 
against finite-size effects. In particular we find that dif¬ 
fusive behaviour also persists for a finite amount of time 
below in the reversible phase, as also described above for 
the one-cycle particle displacement A(r, T). It is how- 
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ever more difficult to measure D in this region, because 
a careful determination of D requires taking the long¬ 
time limit, which is not possible by construction in the 
metastable region. We conclude therefore that the dis¬ 
continuous behaviour of D in Fig.[3fb) appears less con¬ 
vincing than the one of A (r, T) shown in Fig. [2jb) , but 
the overall phenomenology reported in this work appears 
inconsistent with a continuous transition. 


V. NO CHANGE IN MICROSCOPIC 
STRUCTURE 


It was recently argued that the yielding transition in 
oscillatory shear can be detected through the static struc¬ 
ture of the system [27|. Such a behaviour would dif¬ 
fer qualitatively from our conclusion that yielding is re¬ 
vealed through the dynamic evolution of the system. Our 
analysis of the pair correlation function across the yield¬ 
ing transition did not reveal any change in the static 
properties of the system in the two phases, which seems 
to contradict the results of Ref. [23]. To reinforce this 
conclusion, we have measured the exact same quan¬ 
tity that was detected experimentally. In detail, we re¬ 
solve the radial dependence of the static structure factor 
S(q) = jj{p{q)p(-q)) in the (x, z) plane [27], where p(q) 
is the Fourier component of the density at wavevector q. 
Thus we define the angle a from tana = (q z /qx), and fol¬ 
low the a-dependence of the static structure, as proposed 


To obtain statistically reliable data close to the dy¬ 
namic transition, we perform an extensive time average 
over 100 well separated times for the diffusive phase at 
7 o = 0.12. For the non-diffusive phase at 70 = 0.10, time 
averages are not useful and we obtain instead 100 in¬ 
dependent configurations starting from independent ini¬ 
tial conditions. Errorbars are defined from the result¬ 
ing sample-to-sample fluctuations. Because the stress is 
close to a sinusoidal form, we measure the structure either 
when the stress is zero (‘undeformed’ states) and when it 
is maximal a(t) = ±<to (‘deformed’ states). Thus, we ob¬ 
tain 4 measures of the structure at two shear amplitudes, 
for both deformed and undeformed states. 

In Figure |U(a), we show the ^-dependence of the av¬ 
eraged structure factor S(q = |g|) for these 4 cases and 
for N = 1000. We observe that neither the particle re¬ 
versibility nor the deformation seem to affect much the 
structure factor. We now resolve the angular depen¬ 
dence of S(qo,a) using wavevectors in the (x, z) plane 
having an amplitude close to the first peak in the range 
go = 5.9 — 7.0. In Fig. |U(b), we show that the a- 
dependence of the structure factor for the 4 situations de¬ 
fined above is essentially inexistent. Importantly, we do 
not observe any difference for reversible and irreversible 
regimes across yielding. In particular we do not observe 
the oscillations that were detected in the experiments for 
the arrested phase. Furthermore, we also checked that 
average values of other static quantities (such as the en- 



a (rad) 


FIG. 4. Averaged static structure factors for ‘deformed’ 
(when the stress is zero) and undeformed (when the stress 
is maximal) states in the reversible phase at 70 = 0.10 and in 
the diffusive phase at 70 = 0.12, for N = 10 3 . (a) Wavevector 
dependence of the structure with spherical average over all di¬ 
rections. (b) Angular dependence of the structure factor for 
| g| near the first peak as a function of the angle a defined the 
( x,z ) components of q. Errorbars are shown for 70 = 0.10, 
to indicate that our relative resolution is very good (about 2 
%). 


ergy density and pair correlation functions) are similarly 
insensitive to the underlying dynamic transition. These 
conclusions contrast with the results in Ref. [13], which 
are perhaps due to the larger shear rates employed in 
the experiment. Another major contradiction with that 
work is our finding that yielding does not correspond to 
the crossing point of G' and G". Our results show that 
yielding is best interpreted as a loss of reversibility in the 
particle trajectories, which is a purely dynamical con¬ 
cept. 


VI. CONCLUSION 


Together, our results suggest that the yielding transi¬ 
tion of jammed solids under large-amplitude oscillatory 
shear is accompanied by a first-order non-equilibrium 
phase transition, rather than a continuous one. It marks 
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the abrupt emergence of irreversible non-affine particle 
motion. The characteristic strain amplitudes obtained 
from rheology (y x , 7 P i) and from microscopic dynamics 
(bdyn, 7c) are compiled in Fig. [Qc). To a good approxi¬ 
mation, we find that 7 P i « 7d yn ~ 7 C , whereas 7 X is sig¬ 
nificantly larger, corresponding to a large amount of dissi¬ 
pated energy. The dynamic phase transition revealed by 
the discontinuous evolution of single-particle dynamics 
produces a smooth crossover in mechanical properties at 
a critical strain amplitude that appears unrelated to the 
crossing of G"(w) and G"(w). Our conclusions contrast 
with earlier claims of a continuous transition [jj|, UK, 
but appear in very good agreement with observations in 
a sheared emulsion 0. We hope our study will trig¬ 
ger further work in a broader variety of numerical and 


experimental systems to fully establish its generality. 
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